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Abstract: Conceptual engineering design and optimization of laser-based 
imaging techniques and optical diagnostic systems used in the field of 
biomedical optics requires a clear understanding of the light-tissue 
interaction and peculiarities of localization of the detected optical radiation 
within the medium. The description of photon migration within the turbid 
tissue-like media is based on the concept of radiative transfer that forms a 
basis of Monte Carlo (MC) modeling. An opportunity of direct simulation 
of influence of structural variations of biological tissues on the probing light 
makes MC a primary tool for biomedical optics and optical engineering. 
Due to the diversity of optical modalities utilizing different properties of 
light and mechanisms of light-tissue interactions a new MC code is typically 
required to be developed for the particular diagnostic application. In current 
paper introducing an object oriented concept of MC modeling and utilizing 
modern web applications we present the generalized online computational 
tool suitable for the major applications in biophotonics. The computation is 
supported by NVIDEA CUDA Graphics Processing Unit providing 
acceleration of modeling up to 340 times. 
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1. Introduction 

Biophotonics and Biomedical Optics are well established and extremely important research 
areas traditionally associated with biomedical engineering and life sciences. A number of 
achievements in these fields arose due to a successful integration with optical diagnostic (OD) 
and optical/laser based imaging modalities. Nowadays OD technologies are widely used in a 
number of applications including tissue engineering, vascular and developmental biology, 
physiology, morphology, cancer research, material sciences, etc. [1-4]. OD provides a broad 
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variety of practical solutions for clinical diagnostic and therapy in a range of studies from 
single cells and molecules to non-invasive biopsy of specific biological tissues and whole 
organs. The development of OD techniques for specific application is based on the practical 
implementation of laser/optical sources, detectors and various probe configurations. 
Conceptual design, further optimization and/or modification of the particular experimental 
OD systems requires careful selection of various technical parameters, including intensity, 
profile, coherence, polarization, wavelength of incident optical radiation; size, position, 
numerical aperture, sensitivity of the detector; as well as an estimation of how the spatial and 
temporal structural alterations in biological tissues can be linked with and can be 
distinguished by variations of these parameters. Therefore, an accurate realistic description of 
optical radiation propagation within complex biological media is required. 

The description of optical radiation propagation within random media is based on the 
radiative transfer [5] that forms a basis of Monte Carlo (MC) modeling of photons migration 
in biological tissue [6]. Being originally introduced in biomedical optics for the counting of 
fluence rate distribution in biological tissues for the purpose of estimation laser radiation dose 
[7], in the last decades the MC approach has become a primary tool for a number of needs in 
biomedical optics. Incorporated with the computational model of human skin [8] MC 
technique has been used for simulation of skin optical & near-infrared reflectance spectra 
[9,10], analysis of skin fluorescence excitation [11,12], simulation of Optical Coherence 
Tomography (OCT) images of human skin [13,14], analysis of scattering orders and OCT 
image formation [15]. Using a combination of the MC technique and iterative procedure of 
the solution of Bethe-Salpeter equation, the MC approach has been generalized for simulation 
of coherent effects of multiple scattering, such as enhancement of Coherent Back-Scattering 
(CBS) and changes of temporal intensity correlation function depending on the dynamics of 
scattering particles [16,17]. Based on these developments a new approach of handling 
polarization has been introduced and some effects such as a helicity flip of circular 
polarization has been observed [18,19]. The obtained modeling results have been 
comprehensively validated by comparison with the known exact solution by Milne [20,21] 
and with the results of experimental studies of image transfer through the water solution of 
spherical micro-particles of known size and density [22,23]. Meanwhile, a number of other 
MC algorithms has been developed in the past, see for example [24-27]. 

It should be pointed out, however, that due to the diversity of optical modalities utilizing 
different properties of light and mechanisms of light-tissue interactions, including scattering, 
absorption, anisotropy, reflection, refraction, transmittance, interference, polarization, de- 
polarization, coherent scattering, phase shift, time-delaying, etc., a new MC code is typically 
required to be developed for the particular application. In current paper we introduce an object 
oriented concept of MC modeling designed as online computational tool for imitation of 
various aspects of optical/laser radiation propagation within complex disperse multi-layered 
media including biological tissues. 

2. Object oriented concept of Monte Carlo modeling 

Due to a number of practical OD applications the MC model undergoes continuous 
modifications and changes dedicated to inclusion of diverse properties of incident optical/laser 
radiation, configuration of the sources and detectors, structure of the medium and the 
conditions of light detection [8-27]. Past attempts to unify these MC codes are mainly based 
on the use of structured programming [28]. While structured programming is known for years, 
it limits the ability to handle a large code without decreasing its functionality and 
manageability [29]. In practice, the increasing diversity of the MC applications results to a 
substantial growth of the model's source code and leads to the development of a set of 
separate MC codes dedicated each for a particular purpose. 

To generalize and unify the code for a multi-purpose use in various biomedical optics 
applications we apply the Object Oriented Programming (OOP) concept. The OOP is widely 
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used in mainstream application development and has been found extremely effective in design 
of complex multi-parametric systems, providing highly intuitive approach of programming 
[30]. The key features of OOP allow the MC to be separated into logical components, 
described by objects. 

Thus, each photon packets has been defined as the object. The photon objects interact with 
the object - medium or medium components also defined as objects. Splitting the medium into 
the objects allows developing the tissues model more iteratively and uniformly. The 
distribution of scattering centers, macro-inhomogeneities, such as blood vessels, tumors, 
aneurisms, etc. can be formed by combination of 3D elementary volumes (objects) presenting 
spatial variations in the tissues. Moreover, actual structure of a biological tissue can be 
imported into the model as an object (image) provided by OCT, Photo-Acoustic Tomography 
(PAT), ultrasound, MRI, etc. The interaction between object - photons and object - medium 
originates similar as described in [6-23] and omitted for brevity in current paper. 

Utilizing the inheritance feature of OOP the smart hierarchy structure of the code has been 
created to prevent creation of multiple classes for similar tasks. The hierarchy allows 'allied' 
objects to share variables and members, significantly reducing the amount of source code and 
paving the way to extend and generalize the MC for various OD applications. In addition, the 
variations of scattering phase functions, such as Mie, Rayleigh and Henyey-Greenstein [23] 
has been defined using the polymorphism feature of OOP that allows to handle the modeling 
with no changes of the source code. 

Thus, the OOP approach significantly increases the efficiency of the model manageability 
and provides superior opportunities to generalize MC to combine previously developed MC 
models in a way to imitate a particular OD experiment taken into account various features of 
optical radiation and light-tissue interaction. Schematically the structure of the generalized 
MC approach is shown in Fig. 1 . 

As one can see, first, by the selection of a certain application the parameters of the Source, 
Detector and Scattering medium are entered. Depending on the application, the objects are 
tuned to the appropriate feature of light-tissue interaction and the simulation is performed. By 
completing, the output data are prepared in the format utilized in the corresponding OD 
experiment/application. Current MC version supports the applications presented in [6-23] that 
includes estimation of sampling volume [8] offered by a reflectance and/or transmittance 
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Fig. 1. Schematic presentation of the generalized object-oriented MC model structure. Data 
Input and selection of a particular application; Class objects represent exact experimental 
conditions including source-detector geometry & configuration, medium structure, properties 
of incident radiation, etc. 
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geometry (e.g. for fiber optic probe with a small -mm source-detector spacing); simulation of 
reflectance or transmittance optical/near-infrared spectra of the multi-layered media like a 
human skin [9,10]; modeling of skin color depending on volume fraction of melanin and 
blood, blood oxygen saturation; modeling of OCT images with regards coherence and 
polarization properties of probing light (see the details in [13,14]), imitation of spatial 
localization of skin tissues fluorescence excitation [11,12] and simulation of simulation of 
coherent effects of multiple scattering [16,17]. 

Launching of a large number (~10 8 ) of photon packets and computing their interaction 
with medium and with the probe is a highly intensive computational process. Due to the task 
intensity, processing time has always been a significant issue in stochastic modeling, taking a 
few days to complete at the standard CPU. To achieve supreme performance of simulation a 
number programming approaches and optimizations of algorithms have been used in the past, 
including parallel and cluster computing [31,32]. 

We apply recently developed parallel computing framework, known as Compute Unified 
Device Architecture (CUDA), introduced by NVIDEA Corporation. NVIDEA CUDA 
technology provides an unlimited access to computational resources of graphic card: 
processor cores, different types of memory (of various capacity and speed) making Graphics 
Processing Unit (GPU) a massive co-processor in parallel computations [31,32]. 

We utilize the recent, introduced in 2010, CUDA generation, so-called architecture 
codename Fermi. Designed for C/C++ development and easily integrated with the Microsoft 
Visual Studio it makes parallel programming significantly easier, especially in terms of 
project management and debugging. The latest CUDA generation supports most of the OOP 
features like creating classes, inheritance, polymorphism and keeps all dramatic performance 
gains of the CUDA computing. 

The OOP MC model has been developed using CUDA 3.2 C/C++ and supports multiple 
GPUs, providing user's interface and delivering the modeling output. The hardware is 
presented by two GeForce GTX 480 graphic cards with NVIDIA CUDA computing capability 
2.0 totally having up to 2.7 Tflops of computational power on board (single precision). The 
graphic chip is logically divided in hundreds of CUDA cores (in our case 480 per card) 
schematically presented in Fig. 2. Therefore, it executes up to twenty thousands of threads 
simultaneously, without context switch performance losses and very fast (up to 4 Gbit/s) on- 
chip GDDR memory. Graphic chip's shared memory has been used to store the intermediate 
results; constant memory is applied for data input, whereas the global memory is used to store 
parameters of photon objects (e.g. path-length, state of polarization, outlet angles etc.). The 
tiling and cutoff techniques are used to process large data sets and avoid memory bandwidth 
limitations [31,32]. Such a design allows simulating propagation of thousands (15 to 20 




Fig. 2. Schematic presentation of used GPU logical divided on hundreds of independent cores 
allowing creation of thousands of lightweight parallel threads. 
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depending on application and detector parameters) photon objects in the medium 
simultaneously that provides an opportunity for direct imitation of image or wave front 
transfer in the scattering medium [22] or OCT image modeling [13,14]. 

Specially designed for 3D graphic applications mainly for computer games and 
professional software designing this cutting edge graphic technology also incorporates a 
powerful set of instruments applied for optimized simulation of objects motion, rotation, 
reflection, ray-tracing, etc. The NVIDEA CUDA provides GPU-accelerated mathematical 
libraries, such as CULA, CUBLAS — Linear Algebra, CUFFT — Fast Furrier Transform, 
CURAND — Random Number Generators [33]. Their incorporation into MC allows speeding 
up the simulation of each photon packet up to 10 3 times. 

3. Online solution and web integration 

A conceptual design of the online solution is schematically presented on Fig. 3. 



Input / Output 

Load Balancing 

Management 




.NET Interoperability 

COM Integration 
GPU-Web integration 

| Monte Carlo ModellingServer 



Fig. 3. Schematic presentation of the key components of the online MC solution. The server 
hosts a web frontend, which accepts user simulation requests and retrieves results. The 
developed components provide interoperability between the interactive user interface and 
GPUs. 

The MC modeling server provides major hardware acceleration by CUDA supporting 
GPUs. The developed server software consists of the web frontend, management part, GPU- 
Web integration and the CUDA-based MC core described above. 

The frontend of the MC is developed by using Microsoft ASP.NET and Microsoft 
Silverlight technologies [34-37]. Microsoft Silverlight is used for creation of a cross-platform 
lightweight interactive interface to access a particular MC application [37], whereas ASP.NET 
is employed to meet modern design and security requirements [36]. 

The management part of the MC modeling server consists of the Input / Output (I/O) and 
Load Balancing systems. The I/O is created to accept and validate user credentials and online 
input of modeling parameters; to log and store computation enquires and output data into a 
server database; to deliver the results of simulation back to the frontend. The I/O works in a 
tandem with the Load Balancing which allows to manage the server load, monitor running 
simulations and to check the GPUs availability. To make the online MC less vulnerable to the 
common network threats and overloads, the management part has been developed utilizing 
recently added .NET 4.0 features, such as Windows Presentation Foundation (WPF) and 
Windows Communication Foundation (WCF) [38,39]. 

It should be pointed out that .NET solutions are executed in the Common Language 
Runtime (CLR), the core of .NET [35]. Known as a safe managed programming environment 
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the CLR has no direct access to GPUs and other hardware resources. Therefore, to provide the 
interoperability between the web frontend and graphic cards the GPU-Web integration 
component has been developed. The Component Object Model (COM) [40] and .NET 
interoperability have been used to provide interaction between the binaries produced by C# 
and CUDA compilers. The proposed solution allows multiple users to access GPUs 
computational resources online within a reasonable time (in a range from a few seconds to 
several minutes). 

The developed Online Object Oriented MC (0 3 MC) computational tool is now available at 
http://biophotonics.otago.ac.nz through the KAREN network [41], and can be used for various 
applications in biomedical optics as well as in biophotonics and optical engineering. 

4. Results and discussion 

Figure 4-a presents the interactive interface used for the selection of a particular application, 
and Figs. 4-b, 4-c, and 4-d show, respectively, the pages setting the parameters of the medium, 
probe configuration and the outline XYZ boundaries for sampling volume modeling. 

The examples of results of sampling volume and skin reflectance spectra modeling 
counted by 0 3 MC tool for a fiber-optic probe developed in house and used for the 
measurements of optical/near-infrared reflectance spectra of human skin in vivo are presented 
in Fig. 5. Respectively, XZ (depth profile) and XY (surface) images of sampling volume for a 
particular position of source and detector are presented (see Fig. 5). The sampling volume 
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Fig. 4. The frontend of O MC computational tool: (a) interactive user interface providing the 
selection of applications; (b), (c) and (d) are, respectively, the web-pages that allow 
customizing parameters of the medium, fiber-optic probe and the boundaries of sampling 
volume. Available online at: http://biophotonics.otago.ac.nz/MCOnline.aspx. 
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Fig. 5. Example of sampling volume simulation for a fiber-optic probe applied for the 
measurements of human skin reflectance spectra by 0 3 MC: side XZ (left) and surface XY 
(right) projections. Available online at http://biophotonics.otago.ac.nz/MCOnline.aspx. 



presents a spatial localization of the detected signal within the skin tissues providing an 
estimation of which vascular bed is primarily responsible for the measured signal changes. 
The center-to-center source detector separation between fibers has been set as 400 /xm, the 
beam profile has default as Gaussian that is typical for optical fibers. The source and detector 
areas can be overlapped or be coincided. The size and beam profile of the light source, 
numerical aperture, size and position of the detector are taken into account. YZ depth profile 
can be plotted as well. 

Similar interface design and MATLAB style format of the results output are used in other 
applications included modeling human skin color, OCT images, transmittance spectra, CBS, 
confocal probing, etc. The results of simulation of spectra for multi-layered media and skin 
color are similar to the results presented in [9,10]; OCT images modeling and fluorescence are 
based, respectively, on [13,14] and [11,12]. Figure 6 presents an example of the results of 
human skin reflectance spectrum simulation and the results of skin color modeling with 
0 3 MC. The optical parameters of skin are taken same as reported in [9,10]. 

Spectrum of Human Skin 3 Color of Human Skin t3 




Fig. 6. Examples of reflectance spectrum simulation for Type I/II skin, melanin concentration 
2% and normal concentration of blood (left) and skin color counting for Type V/VI skin, 
melanin concentration 25% and normal concentration of blood (right). Available online at 
http://biophotonics.otago.ac.nz/MCOnline.aspx. 

The number of detected photon packets used in the simulation can be varied in a range of 
10 5 -10 8 . Utilizing two GeForce GTX 480 graphic cards the 0 3 MC requires 15.2 seconds to 
simulate the input of 10 8 photon packets. The time of simulation does not depend on 
absorption of the medium, but can be influenced by source-detector geometry increasing from 
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a few seconds to 10 minutes. The results can be saved as an image or downloaded as a data 
file and then used in different scientific graphing software packages. 

The developed 0 3 MC approach has been comprehensively validated through the 
comparison of the results of diffuse reflectance Rd simulation with the known analytical 
results for semi-infinite medium tabulated by Giovanelli [42], with the results of adding- 
doubling method by Prahl [43] and with the results of MCML code developed by Wang et al. 
[7]. The following parameters of the medium have been used in the simulation: fi s = 9 mm~\ 
fi a = 1 mm' 1 , g = 0 (isotropic scattering), n = 1.5. The results of comparison are presented in 
Table 1. Ten simulations with 10 6 photon packets each have been utilized in each case. 
Increasing number of photon packets to 10 8 reduces the error in 0 3 MC for -30%. 



Table 1. Comparison of the diffuse reflectance Rd given by analytical results tabulated by 
Giovanelli [42], the results of the adding-doubling method by Prahl [43], and the results 
of MCML code developed by Wang et al. [7] versus the developed 0 3 MC code 



Reference 


Rd 


Error 


Giovanelli (1955) 


0.2600 




MCML, Wang et al. (1995) 


0.25907 


0.00170 


Adding-doubling, Prahl (1995) 


0.26079 


0.00079 


0 3 MC 


0.25957 


0.00043 



In addition, we compare the results of 0 3 MC simulation of diffuse reflectance Rd and 
transmittance Td for a scattering slab illuminated by normally incident collimated light with 
the results tabulated by van de Hulst [44], with the results of adding-doubling method 
obtained by Prahl [43], and with the results of MCML developed by Wang et al. [7]. The 
following parameters of slab have been used in the simulation: /i s = 9 mm~ l , fi a = 1 mm' 1 , g = 
0.75, n = 1, slab thickness d = 0.2 mm. The results of comparison are listed in Table 2. 



Table 2. Comparison of the diffuse reflectance Rd and transmittance Td given by 
tabulated data by van de Hulst [44], results of the adding-doubling method by Prahl [43], 
and the results of MCML developed by Wang et al. [7] versus the developed 0 3 MC code. 



Reference 


Rd 


error 


Td 


error 


van de Hulst (1957) 


0.09739 




0.66096 




MCML, Wang et al. (1995) 


0.09734 


0.00035 


0.66096 


0.0002 


Adding-doubling, Prahl (1995) 


0.09711 


0.00033 


0.66159 


0.00049 


0 3 MC 


0.09741 


0.00027 


0.66096 


0.00017 



The results presented in Tables 1 and 2 demonstrate a good accuracy of the developed 
code compare to known theoretical predictions and well agreed with the results of MC codes 
developed earlier [7,26,27]. With the further developments of 0 3 MC the fluence rate, 
speckles, intensity/field correlation functions, fluorescence spectra, images of Photo-Acoustic 
Tomography and OCT, pulses propagation, image transfers and other applications currently 
used in biophotonics and biomedical optics will be added to online simulation at 
http://biophotonics.otago.ac.nz/MCOnline.aspx. 
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